1: Ordination

#Load Packages and Prepare Data

library("phyloseq"); packageVersion("phyloseq")

## [1] '1.32.0'

data(GlobalPatterns)
library("ggplot2"); packageVersion("ggplot2")

## [1] '3.3.2'

library("plyr")

theme_set(theme_bw())

Filter low-occurence OTUs

GP = GlobalPatterns
wh0 = genefilter_sample(GP, filterfun_sample(function(x) x > 5), A=0.5*nsamples(GP))
GP1 = prune_taxa(wh0, GP)

Even sampling Depth

GP1 = transform_sample_counts(GP1, function(x) 1E6 * x/sum(x))

Only Keep Five Mos Abundant Phyla

phylum.sum = tapply(taxa_sums(GP1), tax_table(GP1)[, "Phylum"], sum, na.rm=TRUE)
top5phyla = names(sort(phylum.sum, TRUE))[1:5]
GP1 = prune_taxa((tax_table(GP1)[, "Phylum"] %in% top5phyla), GP1)

Create human vs nonhuman categorical variable

human = get_variable(GP1, "SampleType") %in% c("Feces", "Mock", "Skin", "Tongue")
sample_data(GP1)$human <- factor(human)

Four Main ORdination Plots

1) Just OTUs

GP.ord <- ordinate(GP1, "NMDS", "bray")

## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.1333468 
## Run 1 stress 0.1648238 
## Run 2 stress 0.1577391 
## Run 3 stress 0.1333468 
## ... Procrustes: rmse 3.039838e-06  max resid 7.386318e-06 
## ... Similar to previous best
## Run 4 stress 0.161634 
## Run 5 stress 0.2544504 
## Run 6 stress 0.1333468 
## ... Procrustes: rmse 1.045347e-06  max resid 1.726379e-06 
## ... Similar to previous best
## Run 7 stress 0.157644 
## Run 8 stress 0.1748452 
## Run 9 stress 0.153803 
## Run 10 stress 0.2838378 
## Run 11 stress 0.1333468 
## ... Procrustes: rmse 1.722968e-05  max resid 3.041404e-05 
## ... Similar to previous best
## Run 12 stress 0.1622613 
## Run 13 stress 0.1333468 
## ... Procrustes: rmse 2.197773e-06  max resid 6.541958e-06 
## ... Similar to previous best
## Run 14 stress 0.1333468 
## ... New best solution
## ... Procrustes: rmse 4.018371e-06  max resid 1.157531e-05 
## ... Similar to previous best
## Run 15 stress 0.1674589 
## Run 16 stress 0.1333468 
## ... Procrustes: rmse 5.693894e-06  max resid 1.688925e-05 
## ... Similar to previous best
## Run 17 stress 0.1333468 
## ... Procrustes: rmse 3.422433e-06  max resid 9.37327e-06 
## ... Similar to previous best
## Run 18 stress 0.146027 
## Run 19 stress 0.1529609 
## Run 20 stress 0.1625845 
## *** Solution reached

p1 = plot_ordination(GP1, GP.ord, type="taxa", color="Phylum", title="taxa")
print(p1)

Simplify Visually

p1 + facet_wrap(~Phylum, 3)

2) Just Samples

p2 = plot_ordination(GP1, GP.ord, type="samples", color="SampleType", shape="human") 
p2 + geom_polygon(aes(fill=SampleType)) + geom_point(size=5) + ggtitle("samples")

3) Biplot Graphic

p3 = plot_ordination(GP1, GP.ord, type="biplot", color="SampleType", shape="Phylum", title="biplot")
# Some stuff to modify the automatic shape scale
GP1.shape.names = get_taxa_unique(GP1, "Phylum")
GP1.shape <- 15:(15 + length(GP1.shape.names) - 1)
names(GP1.shape) <- GP1.shape.names
GP1.shape["samples"] <- 16
p3 + scale_shape_manual(values=GP1.shape)

4) Split Graphic

p4 = plot_ordination(GP1, GP.ord, type="split", color="Phylum", shape="human", label="SampleType", title="split") 
p4

Set to Default Color Scale

gg_color_hue <- function(n){
    hues = seq(15, 375, length=n+1)
    hcl(h=hues, l=65, c=100)[1:n]
}
color.names <- levels(p4$data$Phylum)
p4cols <- gg_color_hue(length(color.names))
names(p4cols) <- color.names
p4cols["samples"] <- "black"
p4 + scale_color_manual(values=p4cols)

Supported Ordination Methods

dist = "bray"
ord_meths = c("DCA", "CCA", "RDA", "DPCoA", "NMDS", "MDS", "PCoA")
plist = lapply(as.list(ord_meths), function(i, physeq, dist){
        ordi = ordinate(physeq, method=i, distance=dist)
        plot_ordination(physeq, ordi, "samples", color="SampleType")
}, GP1, dist)

## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.1333468 
## Run 1 stress 0.1385323 
## Run 2 stress 0.156695 
## Run 3 stress 0.1333468 
## ... New best solution
## ... Procrustes: rmse 1.296711e-06  max resid 3.456084e-06 
## ... Similar to previous best
## Run 4 stress 0.1617642 
## Run 5 stress 0.1626327 
## Run 6 stress 0.1333468 
## ... Procrustes: rmse 1.762521e-06  max resid 4.140732e-06 
## ... Similar to previous best
## Run 7 stress 0.1522771 
## Run 8 stress 0.1529686 
## Run 9 stress 0.146027 
## Run 10 stress 0.1820691 
## Run 11 stress 0.1569043 
## Run 12 stress 0.1625845 
## Run 13 stress 0.146027 
## Run 14 stress 0.1333468 
## ... Procrustes: rmse 3.472052e-06  max resid 1.03564e-05 
## ... Similar to previous best
## Run 15 stress 0.1529686 
## Run 16 stress 0.3021346 
## Run 17 stress 0.1449078 
## Run 18 stress 0.17179 
## Run 19 stress 0.1333468 
## ... Procrustes: rmse 8.106219e-06  max resid 2.431961e-05 
## ... Similar to previous best
## Run 20 stress 0.1333468 
## ... New best solution
## ... Procrustes: rmse 1.75444e-06  max resid 4.440386e-06 
## ... Similar to previous best
## *** Solution reached

names(plist) <- ord_meths

pdataframe = ldply(plist, function(x){
    df = x$data[, 1:2]
    colnames(df) = c("Axis_1", "Axis_2")
    return(cbind(df, x$data))
})
names(pdataframe)[1] = "method"

p = ggplot(pdataframe, aes(Axis_1, Axis_2, color=SampleType, shape=human, fill=SampleType))
p = p + geom_point(size=4) + geom_polygon()
p = p + facet_wrap(~method, scales="free")
p = p + scale_fill_brewer(type="qual", palette="Set1")
p = p + scale_colour_brewer(type="qual", palette="Set1")
p

plist[[2]]

p = plist[[2]] + scale_colour_brewer(type="qual", palette="Set1")
p = p + scale_fill_brewer(type="qual", palette="Set1")
p = p + geom_point(size=5) + geom_polygon(aes(fill=SampleType))
p

MDS (“PCoA”) on Unifrac Distances

ordu = ordinate(GP1, "PCoA", "unifrac", weighted=TRUE)
plot_ordination(GP1, ordu, color="SampleType", shape="human")

p = plot_ordination(GP1, ordu, color="SampleType", shape="human")
p = p + geom_point(size=7, alpha=0.75)
p = p + scale_colour_brewer(type="qual", palette="Set1")
p + ggtitle("MDS/PCoA on weighted-UniFrac distance, GlobalPatterns")

2: Alpha Diversity Graphics

Load packages, set parameters

library("phyloseq"); packageVersion("phyloseq")

## [1] '1.32.0'

data("GlobalPatterns")

library("ggplot2"); packageVersion("ggplot2")

## [1] '3.3.2'

theme_set(theme_bw())
pal = "Set1"
scale_colour_discrete <-  function(palname=pal, ...){
  scale_colour_brewer(palette=palname, ...)
}
scale_fill_discrete <-  function(palname=pal, ...){
  scale_fill_brewer(palette=palname, ...)
}

Prepare Data

GP <- prune_species(speciesSums(GlobalPatterns) > 0, GlobalPatterns)

## Warning: 'prune_species' is deprecated.
## Use 'prune_taxa' instead.
## See help("Deprecated") and help("phyloseq-deprecated").

## Warning: 'speciesSums' is deprecated.
## Use 'taxa_sums' instead.
## See help("Deprecated") and help("phyloseq-deprecated").

Plot Examples

plot_richness(GP)

plot_richness(GP, measures=c("Chao1", "Shannon"))

plot_richness(GP, x="SampleType", measures=c("Chao1", "Shannon"))

sampleData(GP)$human <- getVariable(GP, "SampleType") %in% c("Feces", "Mock", "Skin", "Tongue")

## Warning: 'getVariable' is deprecated.
## Use 'get_variable' instead.
## See help("Deprecated") and help("phyloseq-deprecated").

## Warning: 'sampleData' is deprecated.
## Use 'sample_data' instead.
## See help("Deprecated") and help("phyloseq-deprecated").

## Warning: 'sampleData<-' is deprecated.
## Use 'sample_data<-' instead.
## See help("Deprecated") and help("phyloseq-deprecated").

plot_richness(GP, x="human", color="SampleType", measures=c("Chao1", "Shannon"))

GPst = merge_samples(GP, "SampleType")
# repair variables that were damaged during merge (coerced to numeric)
sample_data(GPst)$SampleType <- factor(sample_names(GPst))
sample_data(GPst)$human <- as.logical(sample_data(GPst)$human)

p = plot_richness(GPst, x="human", color="SampleType", measures=c("Chao1", "Shannon"))
p + geom_point(size=5, alpha=0.7)

More about GGPLOT2

Check the lists present

p$layers

## [[1]]
## geom_point: na.rm = TRUE
## stat_identity: na.rm = TRUE
## position_identity 
## 
## [[2]]
## mapping: ymax = ~value + se, ymin = ~value - se 
## geom_errorbar: na.rm = FALSE, orientation = NA, width = 0.1, width = 0.1, flipped_aes = FALSE
## stat_identity: na.rm = FALSE
## position_identity

p$layers <- p$layers[-1]
p + geom_point(size=5, alpha=0.7)

3: Heatmap Plots

Examples

library("phyloseq"); packageVersion("phyloseq")

## [1] '1.32.0'

library("ggplot2"); packageVersion("ggplot2")

## [1] '3.3.2'

theme_set(theme_bw())

Plot a 300-Taxa Dataset

data("GlobalPatterns")
gpt <- subset_taxa(GlobalPatterns, Kingdom=="Bacteria")
gpt <- prune_taxa(names(sort(taxa_sums(gpt),TRUE)[1:300]), gpt)
plot_heatmap(gpt, sample.label="SampleType")

## Warning: Transformation introduced infinite values in discrete y-axis

Subset a smaller dataset based on an Archaeal phylum

gpac <- subset_taxa(GlobalPatterns, Phylum=="Crenarchaeota")

Default plot_heatmap settings

plot_heatmap(gpac)

## Warning: Transformation introduced infinite values in discrete y-axis

Re-label by a sample variable and taxonomic family

(p <- plot_heatmap(gpac, "NMDS", "bray", "SampleType", "Family"))

## Warning: Transformation introduced infinite values in discrete y-axis

Relabel Axis Titles

p$scales$scales[[1]]$name <- "My X-Axis"
p$scales$scales[[2]]$name <- "My Y-Axis"
print(p)

## Warning: Transformation introduced infinite values in discrete y-axis

Now repeat the plot, but change the color scheme

plot_heatmap(gpac, "NMDS", "bray", "SampleType", "Family", low="#000033", high="#CCFF66")

## Warning: Transformation introduced infinite values in discrete y-axis

plot_heatmap(gpac, "NMDS", "bray", "SampleType", "Family", low="#000033", high="#FF3300")

## Warning: Transformation introduced infinite values in discrete y-axis

plot_heatmap(gpac, "NMDS", "bray", "SampleType", "Family", low="#000033", high="#66CCFF")

## Warning: Transformation introduced infinite values in discrete y-axis

Dark on White Color Scheme

plot_heatmap(gpac, "NMDS", "bray", "SampleType", "Family", low="#66CCFF", high="#000033", na.value="white")

## Warning: Transformation introduced infinite values in discrete y-axis

Cream Colored

plot_heatmap(gpac, "NMDS", "bray", "SampleType", "Family", low="#FFFFCC", high="#000033", na.value="white")

## Warning: Transformation introduced infinite values in discrete y-axis

Different Ordination Methods, Distances

Jaccard Distance

plot_heatmap(gpac, "NMDS", "jaccard")

## Warning: Transformation introduced infinite values in discrete y-axis

Detrended Correspodence Analysis

plot_heatmap(gpac, "DCA", "none", "SampleType", "Family")

## Warning: Transformation introduced infinite values in discrete y-axis

Unconstrained Reduncancy Analysis

plot_heatmap(gpac, "RDA", "none", "SampleType", "Family")

## Warning: Transformation introduced infinite values in discrete y-axis

PCoA/MDS of Bray-Curtis Distance

plot_heatmap(gpac, "PCoA", "bray", "SampleType", "Family")

## Warning: Transformation introduced infinite values in discrete y-axis

PCoA/MDS of Unifrac Distance

plot_heatmap(gpac, "PCoA", "unifrac", "SampleType", "Family")

## Warning: Transformation introduced infinite values in discrete y-axis

Weighted-UniFrac distance and MDS/PCoA ordination

plot_heatmap(gpac, "MDS", "unifrac", "SampleType", "Family", weighted=TRUE)

## Warning: Transformation introduced infinite values in discrete y-axis

Hierarchichal Clustering

heatmap(otu_table(gpac))

4: Networks

Load Package and Example Data

library(phyloseq); packageVersion("phyloseq")

## [1] '1.32.0'

packageVersion("ggplot2")

## [1] '3.3.2'

data(enterotype)

set.seed(711L)

enterotype = subset_samples(enterotype, !is.na(Enterotype))

The Plot_Net Function

plot_net(enterotype, maxdist = 0.4, point_label = "Sample_ID")

Map as Color and Shape

plot_net(enterotype, maxdist = 0.3, color = "SeqTech", shape="Enterotype")

The Plot_Network Function

ig <- make_network(enterotype, max.dist=0.3)

plot_network(ig, enterotype)

Map Samples Variables as Color and Shape

plot_network(ig, enterotype, color="SeqTech", shape="Enterotype", line_weight=0.4, label=NULL)

Lower Maximum Distance

ig <- make_network(enterotype, max.dist=0.2)
plot_network(ig, enterotype, color="SeqTech", shape="Enterotype", line_weight=0.4, label=NULL)

Replace Jaccard with Bray-Curtis

ig <- make_network(enterotype, dist.fun="bray", max.dist=0.3)
plot_network(ig, enterotype, color="SeqTech", shape="Enterotype", line_weight=0.4, label=NULL)

Lab11x CKitrinos Markdown

1: Ordination

Filter low-occurence OTUs

Even sampling Depth

Only Keep Five Mos Abundant Phyla

Create human vs nonhuman categorical variable

Four Main ORdination Plots

1) Just OTUs

Simplify Visually

2) Just Samples

3) Biplot Graphic

4) Split Graphic

Set to Default Color Scale

Supported Ordination Methods

MDS (“PCoA”) on Unifrac Distances

2: Alpha Diversity Graphics

Load packages, set parameters

Prepare Data

Plot Examples

More about GGPLOT2

Check the lists present

3: Heatmap Plots

Examples

Plot a 300-Taxa Dataset

Subset a smaller dataset based on an Archaeal phylum

Default plot_heatmap settings

Re-label by a sample variable and taxonomic family

Relabel Axis Titles

Now repeat the plot, but change the color scheme

Dark on White Color Scheme

Cream Colored

Different Ordination Methods, Distances

Jaccard Distance

Detrended Correspodence Analysis

Unconstrained Reduncancy Analysis

PCoA/MDS of Bray-Curtis Distance

PCoA/MDS of Unifrac Distance

Weighted-UniFrac distance and MDS/PCoA ordination

Hierarchichal Clustering

4: Networks

Load Package and Example Data

The Plot_Net Function

Map as Color and Shape

The Plot_Network Function

Map Samples Variables as Color and Shape

Lower Maximum Distance

Replace Jaccard with Bray-Curtis